System And Method For Comparing Dental X-Ray Images

ABSTRACT

A system for comparing dental X-ray images includes a positional displacement calculator calculating a positional displacement between dental X-ray test and reference images by using phase-only correlation, a positional displacement corrector correcting the positional displacement, a base point extractor defining, as a base image, any one of the dental X-ray test and reference images, and defining, as a corresponding image, the other one of the two dental images, and extracting base points from the base image, a corresponding point extractor extracting corresponding points, which correspond to the base points, from the corresponding image, a correspondence calculator calculating correspondence between the base points and the corresponding points, a nonlinear distortion corrector correcting a nonlinear distortion between the base image and the corresponding image, based on the correspondence, and a similarity calculator finding, by using phase-only correlation, a similarity between the base image and the corresponding image.

CROSS REFERENCE TO RELATED APPLICATIONS AND INCORPORATION BY REFERENCE

This application is based upon and claims the benefit of priority fromprior Japanese Patent Application P2008-120622 filed on May 2, 2008; theentire contents of which are incorporated by reference herein.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a method for comparing images, andparticularly to a system and method for comparing dental X-ray images.

2. Description of the Related Art

Dental radiographs are used for the purpose of finding a slight changein the internal structure of a bone, monitoring the progress of adisease, developing treatment policies, or determining an identity (seeLehmann, T. M., Grondahl, H. G., & Benn, D. K. (2000). Computer-basedregistration for digital subtraction in dental radiology.Dentomaxillofacial Radiology, 29, 323-346., for example). When dentalradiographs are used for any one of these purposes, a radiograph takenseveral weeks ago or several years ago needs to be accurately comparedwith a radiograph taken recently. In order to realize a computer-aideddiagnosis (CAD) in a dental treatment, it is important that a dentalradiograph taken recently should be accurately positioned to a dentalradiograph taken in the past. However, because dental radiographs aretaken with the films placed in the complicated and narrow oral cavity, anonlinear distortion exists between the dental radiograph taken recentlyand the dental radiograph taken in the past, even though the dentalradiographs are taken of the same oral cavity of the same person. Thisbrings about a problem that it is difficult to compare the dentalradiograph taken recently and the dental radiograph taken in the past.

SUMMARY OF THE INVENTION

An aspect of present invention inheres in a system for comparing dentalX-ray images including a positional displacement calculator configuredto calculate a positional displacement between a dental X-ray test imageand a dental X-ray reference image by using phase-only correlation, apositional displacement corrector configured to correct the positionaldisplacement, a base point extractor configured to define any one of thedental X-ray test image whose positional displacement is corrected andthe dental X-ray reference image as a base image, and to define theother one of the two dental images as a corresponding image, as well asto extract multiple base points from the base image, a correspondingpoint extractor configured to extract multiple corresponding points,which correspond to the multiple base points, from the correspondingimage, a correspondence calculator configured to calculatecorrespondence between the multiple base points and the multiplecorresponding points nonlinearly distorted from the multiple basepoints, a nonlinear distortion corrector configured to correct anonlinear distortion between the base image and the corresponding image,based on the correspondence, and a similarity calculator configured tofind a similarity between the base image whose nonlinear distortion iscorrected and the corresponding image by using phase-only correlation.

Another aspect of the present invention inheres in a method forcomparing dental X-ray images including calculating a positionaldisplacement between a dental X-ray test image and a dental X-rayreference image by using phase-only correlation, correcting thepositional displacement, defining any one of the dental X-ray test imagewhose positional displacement is corrected and the dental X-rayreference image as a base image, and defining the other one of the twodental images as a corresponding image, as well as extracting multiplebase points from the base image, extracting multiple correspondingpoints, which correspond to the multiple base points, from thecorresponding image, calculating correspondence between the multiplebase points and the multiple corresponding points nonlinearly distortedfrom coordinates of the multiple base points, correcting a nonlineardistortion between the base image and the corresponding image, based onthe correspondence, and finding a similarity between the base image andthe corresponding image between which the nonlinear distortion iscorrected by using phase-only correlation.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of a system for comparing dental X-rayimages according to an embodiment of the present invention.

FIG. 2 is a schematic diagram showing how to take dental X-ray imagesaccording to the embodiment of the present invention.

FIG. 3 is a first schematic diagram of a tooth which a dental X-rayimage according to the embodiment of the present invention is taken of.

FIG. 4 is a second schematic diagram of the tooth which another dentalX-ray image according to the embodiment of the present invention istaken of.

FIG. 5 is a third schematic diagram of the tooth which yet anotherdental X-ray image according to the embodiment of the present inventionis taken of.

FIG. 6 shows a dental X-ray test image and a dental X-ray referenceimage according to the embodiment of the present invention.

FIG. 7 shows a contrast-corrected dental X-ray test image and acontrast-corrected dental X-ray reference image according to theembodiment of the present invention.

FIG. 8 shows a dental X-ray test image, and a dental X-ray referenceimage, which result from correcting a rotation angle, a scaling ratioand the amount of parallel translation according to the embodiment ofthe present invention.

FIG. 9 shows a base image from which multiple base points are extractedand a corresponding image according to the embodiment of the presentinvention.

FIG. 10 shows the base image and the corresponding image in whichmultiple corresponding points are found according to the embodiment ofthe present invention.

FIG. 11 shows grids on a corresponding image which results fromcorrecting a nonlinear distortion from the base image according to theembodiment of the present invention.

FIG. 12 shows the corresponding image which results from correcting thenonlinear distortion from the base image according to the embodiment ofthe present invention.

FIG. 13 shows a base image and a corresponding image from which commonareas are extracted according to the embodiment of the presentinvention.

FIG. 14 is a flowchart showing a method for comparing dental X-rayimages according to the embodiment of the present invention.

FIG. 15 shows first dental radiographs respectively taken before andafter a dental treatment according to an example of the embodiment ofthe present invention.

FIG. 16 shows second dental radiographs respectively taken before andafter a dental treatment according to the example of the embodiment ofthe present invention.

FIG. 17 shows third dental radiographs respectively taken before andafter a dental treatment according to the example of the embodiment ofthe present invention.

FIG. 18 is a graph showing a cumulative match curve according to theexample of the embodiment of the present invention.

FIG. 19 is a graph showing frequencies of comparison scores according tothe example of the embodiment of the present invention.

FIG. 20 shows a first comparison result according to the example of theembodiment of the present invention.

FIG. 21 shows a second comparison result according to the example of theembodiment of the present invention.

FIG. 22 shows a third comparison result according to the example of theembodiment of the present invention.

FIG. 23 shows a fourth comparison result according to the example of theembodiment of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

Descriptions will be provided hereinbelow for an embodiment of thepresent invention. In descriptions to the following drawings, the sameor similar components are denoted by the same or similar referencenumerals. It should be noted that the drawings are schematic. For thisreason, specific dimensions and the like should be understood with thefollowing descriptions taken into consideration. It goes without sayingthat the dimensional relationships and ratios among some components aredifferent from one drawing to another.

With reference to FIG. 1, a system for comparing dental X-ray imagesaccording to the embodiment of the present invention includes, an imageacquiring device 200, such as a scanner, for acquiring a dental X-raytest image that has been taken, and a central processing unit (CPU) 300connected to the image acquiring device 200. The CPU 300 includes apositional displacement calculator 301 configured to calculate apositional displacement between the dental X-ray test image and a dentalX-ray reference image, based on phase-only correlation, and a positionaldisplacement corrector 302 configured to correct the positionaldisplacement between the dental X-ray test image and the dental X-rayreference image.

In addition, the CPU 300 includes a definer 303 configured to defineanyone of the dental X-ray test image whose positional displacement iscorrected and the dental X-ray reference image as a base image, and todefine the other one of the two dental images as a corresponding image,a base point extractor 304 configured to extract multiple base points ateach of which a pixel value changes beyond a threshold value, from thebase image, and a corresponding point extractor 305 configured toextract multiple corresponding points, which correspond to the multiplebase points, from the corresponding image. A same object is recorded atthe base point and the corresponding point.

Moreover, the CPU 300 includes a correspondence calculator 306configured to calculate correspondence between coordinates of themultiple base points and coordinates of the corresponding points whichare nonlinearly distorted from the coordinates of the multiple basepoints, a nonlinear distortion corrector 307 configured to correct anonlinear distortion between the base image and the corresponding image,based on the correspondence, and a similarity calculator 308 configuredto find a similarity between the base image and the corresponding imagebetween which the nonlinear distortion is corrected, based on thephase-only correlation.

Hereinbelow, consider the dental X-ray test image and the dental X-rayreference image as images each with N₁×N₂ pixels. In addition, thedental X-ray test image is represented with g(n₁, n₂), and the dentalX-ray reference image is represented with f(n₁, n₂).

As shown in FIG. 2, the dental X-ray test image g(n₁, n₂) is taken byplacing a small film 50 behind a tooth 10, and subsequently by placingan X-ray irradiator 100 near the face with the film 50 being held by afinger, as well as thereafter by irradiating an X-ray on the film 50from the X-ray irradiator 100. A previously acquired dental X-rayreference image f(n₁, n₂) is also taken in the same manner as the dentalX-ray test image g(n₁, n₂). For the purpose of taking an idealisticdental X-ray image, it is necessary that the film 50 should be placed,in parallel with the tooth 10, and perpendicularly to the X-rayirradiated on the film 50. However, because the film 50 and the X-rayirradiator 100 are manually placed by a picture taker each time thedental X-ray image is taken, the placement position and angle of each ofthe film 50 and the X-ray irradiator may vary each time the dental X-rayimage is taken. As a result, even though the dental X-ray test imageg(n₁, n₂) and the dental X-ray reference image f(n₁, n₂) were taken ofthe same oral cavity area of the same person, the tooth 10 as the objectin the dental X-ray test image may be parallel translated and rotated incomparison with the tooth 10 as the object in the dental X-ray referenceimage because the dental X-ray test image and the dental X-ray referenceimage were taken at different times.

For example, the tooth 10 in the dental X-ray image, shown in FIG. 4,which is obtained when the X-ray irradiator 100 is placed in a secondposition shown in FIG. 2 may be shorter than the tooth 10 in the dentalX-ray image, shown in FIG. 3, which is obtained when the X-rayirradiator 100 is placed in a first position shown in FIG. 2. Inaddition, the tooth 10 in the dental X-ray image, shown in FIG. 5, whichis obtained when the X-ray irradiator 100 is placed in a third positionshown in FIG. 2 may be far shorter than the tooth 10 in the dental X-rayimage, shown in FIG. 4, which is obtained when the X-ray irradiator 100is placed in the second position shown in FIG. 2.

Furthermore, the distance between the film 50 and the X-ray irradiator100 may vary each time the dental X-ray image is taken of the tooth 10.For this reason, the tooth 10 as the object in the dental X-ray testimage g(n₁, n₂) may be larger or smaller than the tooth 10 as the objectin the dental X-ray reference image f(n₁, n₂), because the dental X-raytest image and the dental X-ray reference image were taken at differenttimes. Moreover, because the film 50 is pressed against the tooth 10with a finger while the dental X-ray image is being taken of the tooth10, the film 50 may be nonlinearly distorted. The nonlinear distortionof the film 50 may vary each time the dental X-ray image is taken.Additionally, the dental X-ray test image g(n₁, n₂) and the dental X-rayreference image f(n₁, n₂) may be blurred because of noise and the likeincluded while the dental X-ray images are being taken. Examples of thedental X-ray test image g(n₁, n₂) and the dental X-ray reference imagef(n₁, n₂) are shown in FIG. 6.

The CPU 300 shown in FIG. 1 further includes a contrast corrector 401.The contrast corrector 401 enhances the contrasts respectively of thedental X-ray test image g(n₁, n₂) and the dental X-ray reference imagef(n₁, n₂), and thus generates the contrast-corrected dental X-ray testimage g_(e)(n₁, n₂) and the contrast-corrected dental X-ray referenceimage f_(e)(n₁, n₂) which are shown in FIG. 7. The contrast corrector401 shown in FIG. 1 is designed to enhance a contrast by using LACE(Local Area Contrast Enhancement) and morphological filter. The contrastenhancement of the dental X-ray test image g(n₁, n₂) and the dentalX-ray reference image f(n₁, n₂) by the contrast corrector 401 makes itpossible to obviously increase the accuracy with which thebelow-described comparison score is calculated.

The positional displacement calculator 301 calculates thetwo-dimensional discrete Fourier transform G_(e)(k₁, k₂) of thecontrast-corrected dental X-ray test image g_(e)(n₁, n₂) by usingEquation (1). In the following example, we assume that the index rangesin the discrete space are n₁=−M₁, . . . ,M₁(M₁>0) and n₂=−M₂, . . . ,M₂(M₂>0) for mathematical simplicity, and hence N₁=2M₁+1 and N₂=2M₂+1.Note that, in the case of the present embodiment, the index ranges inthe discrete space are symmetrical between positive and negativenumbers. And, the number N₁ of vertically-arrayed pixels and the numberN₂ of horizontally-arrayed pixels are odd numbers, for explanatoryconvenience. However, it goes without saying that the two-dimensionaldiscrete Fourier transform may be generalized in order that non-negativeindex ranges can be used, and concurrently in order that the number N₁of vertically-arrayed pixels and the number N₂ of horizontally-arrayedpixels can be respectively set at arbitrary positive integers.G_(e)(k₁,k₂) is given by

$\begin{matrix}\begin{matrix}{{G_{e}\left( {k_{1},k_{2}} \right)} = {\sum\limits_{{n\; 1} = {{- M}\; 1}}^{M\; 1}{\sum\limits_{{n\; 2} = {{- M}\; 2}}^{M\; 2}{{g_{e}\left( {n_{1},n_{2}} \right)}W_{N\; 1}^{k\; 1n\; 1}W_{N\; 2}^{k\; 2n\; 2}}}}} \\{= {{A_{G}\left( {k_{1},k_{2}} \right)}^{j\theta}G^{({{k\; 1},{k\; 2}})}}}\end{matrix} & {{Equation}\mspace{20mu} (1)}\end{matrix}$

where W_(N1)=e^(−j(2π/N1)) and W_(N2)=e^(−j(2π/N2)), and A_(G)(k₁, k₂)is an amplitude component of the dental X-ray test image g_(e)(n₁, n₂),as well as θ_(G)(k₁, k₂) is a phase component of the dental X-ray testimage g_(e)(n₁, n₂).

In addition, the positional displacement calculator 301 calculates thetwo-dimensional discrete Fourier transform F_(e)(k₁, k₂) of thecontrast-corrected dental X-ray reference image f_(e)(n₁, n₂) by use of

$\begin{matrix}\begin{matrix}{{F_{e}\left( {k_{1},k_{2}} \right)} = {\sum\limits_{{n\; 1} = {{- M}\; 1}}^{M\; 1}{\sum\limits_{{n\; 2} = {{- M}\; 2}}^{M\; 2}{{f_{e}\left( {n_{1},n_{2}} \right)}W_{N\; 1}^{k\; 1\; n\; 1}W_{N\; 2}^{k\; 2\; n\; 2}}}}} \\{= {{A_{F}\left( {k_{1},k_{2}} \right)}^{j\theta}F^{({{k\; 1},{k\; 2}})}}}\end{matrix} & {{Equation}\mspace{20mu} (2)}\end{matrix}$

where A_(F)(k₁, k₂) is an amplitude component of the dental X-rayreference image f_(e)(n₁, n₂), and θ_(F)(k₁, k₂) is a phase component ofthe dental X-ray reference image f_(e)(n₁, n₂).

Thereafter, the positional displacement calculator 301 calculates thesquare root of the amplitude spectrum of the dental X-ray test image(|G_(e)(k₁, k₂)|)^(1/2) and the square root of the amplitude spectrum ofthe dental X-ray reference image (|F_(e)(k₁, k₂)|)^(1/2). The amplitudespectra are free from the influence of the amount of paralleltranslation of the object between the dental X-ray test image and thedental X-ray reference image, because the amplitude spectra areinfluenced by nothing but the rotation angle θ and the scaling ratio κbetween the dental X-ray test image and the dental X-ray referenceimage.

In the case of natural images, almost all of the energy concentrates ona lower frequency range. For this reason, information existing in ahigher frequency range generally is enhanced and thus visualized byusing the logarithms of the respective amplitude spectra log {|G_(e)(k₁,k₂)|+1} and log {|F_(e)(k₁, k₂)|+1}. However, in the case of the dentalX-ray images, the amplitude components including information on themagnification, reduction and rotation concentrates on a far lowerfrequency range than in the case of the general natural images. For thisreason, if the amplitude spectra are enhanced by using these logarithms,information with a lower signal-to-noise ratio (SN ratio) is farenhanced. This reduces the accuracy with which the rotation angle θ andthe scaling ratio κ are estimated. By contrast, use of the square rootsof the respective amplitude spectra makes it possible to smoothlyenhance the amplitude spectra of the dental X-ray images from the lowerfrequency range to the higher frequency range.

Subsequently, the positional displacement calculator 301 computes thelog-polar transform (or the Fourier-Mellin transform) of the square root(|G_(e)(k₁, k₂)|)^(1/2) of the amplitude spectrum of the dental X-raytest image, and thus calculates G_(LF)(k₁′, k₂′). In addition, thepositional displacement calculator 301 computes the log-polar transform(or the Fourier-Mellin transform) of the square root (|F_(e)(k₁,k₂)|)^(1/2) of the amplitude spectrum of the dental X-ray test image,and thus calculates F_(LP)(k₁′, k₂′) Where k₁′ and k₂′ are indices ofthe dental X-ray images to which the log-polar transform is applied.

Thereafter, the positional displacement calculator 301 calculates thenormalized cross power spectrum R_(FlpGlp)(k₁′, k₂′) of F_(LP)(k₁′, k₂′)and G_(LP)(k₁′, k₂′). RF_(FlpGlp)(k₁′, k₂′) is given by

$\begin{matrix}\begin{matrix}{{R_{FlpGlp}\left( {k_{1}^{\prime},k_{2}^{\prime}} \right)} = \frac{{F_{Lp}\left( {k_{1}^{\prime},k_{2}^{\prime}} \right)}\overset{\_}{G_{LP}\left( {k_{1}^{\prime},k_{2}^{\prime}} \right)}}{{{F_{LP}\left( {k_{1}^{\prime},k_{2}^{\prime}} \right)}\overset{\_}{G_{LP}\left( {k_{1}^{\prime},k_{2}^{\prime}} \right)}}}} \\{= ^{{j\theta}{({k_{1}^{\prime},k_{2}^{\prime}})}}}\end{matrix} & {{Equation}\mspace{20mu} (3)}\end{matrix}$

where θ(k₁′, k₂′)=θ_(Flp)(k₁′, k₂′) −θ_(Glp)(k₁′, k₂′), andG_(Lp)(k₁′,k₂′) is the complex conjugate of G_(LP)(k₁′, k₂′).

Subsequently, the positional displacement calculator 301 calculates theband-limited phase-only correlation (BLPOC) function r_(FlpGlp)^(K1K2)(n₁, n₂) as the two-dimensional inverse discrete Fouriertransform of the normalized cross power spectrum R_(FlpGlp)(k₁′, k₂′).r_(FlpGlp) ^(K1K2)(n₁, n₂) is given by

$\begin{matrix}{{r_{FlpGlp}^{K\; 1K\; 2}\left( {n_{1},n_{2}} \right)} = {\frac{1}{L_{1}L_{2}}{\sum\limits_{{k\; 1^{\prime}} = {{- K}\; 1}}^{K\; 1}{\sum\limits_{{k\; 2^{\prime}} = {{- K}\; 2}}^{K\; 2}{{R_{FlpGlp}\left( {k_{1}^{\prime},k_{2}^{\prime}} \right)} \times W_{L_{1}}^{{- k_{1}^{\prime}}n_{1}}W_{L_{2}}^{{- k_{2}^{\prime}}n_{2}}}}}}} & {{Equation}\mspace{20mu} (4)}\end{matrix}$

where K₁(0<K₁≦M₁) and K₂(0<K₂≦M₂) are effective bands of thetwo-dimensional inverse discrete Fourier transform, as well as L₁=2K₁+1and L₂=2K₂+1. In addition, for example, K₁/M₁₌0.2 and K₂/M₂=0.2.

The use of the band-limited phase-only correlation function limits thetwo-dimensional inverse discrete Fourier transform of the normalizedcross power spectrum R_(FlpGlp)(k₁′, k₂′) in the effective bands of theimage texture, and thus concentrates the energy of the correlation peak.This makes it possible to enhance the image identification performancewhile eliminating the influence of the high-frequency components with alower reliability. In addition, the use of the band-limited phase-onlycorrelation function makes the size of the two dimensional inversediscrete Fourier transform smaller than the use of the phase-onlycorrelation function, thus reducing the amount of calculation. In spiteof this size reduction, the accuracy with which the positionaldisplacement is estimated by using the band-limited phase-onlycorrelation function is almost equal to the accuracy with which thepositional displacement is estimated by using the phase-only correlationfunction.

The positional displacement calculator 301 detects the location of thecorrelation peak of the band-limited phase-only correlation functionr_(FlpGlp) ^(K1K2)(n₁, n₂), and thus calculates the rotation angle θ ofthe dental X-ray test image g_(e)(n₁, n₂) to the dental X-ray referenceimage f_(e)(n₁, n₂) and the scaling ratio κ of the dental X-ray testimage g_(e)(n₁, n₂) to the dental X-ray reference image f_(e)(n₁, n₂).The positional displacement corrector 302 corrects the dental X-ray testimage g_(e)(n₁, n₂) in order that the rotation angle θ of the dentalX-ray test image g_(e)(n₁, n₂) to the dental X-ray reference imagef_(e)(n₁, n₂) can be equal to 0 (zero), and concurrently in order thatthe scaling ratio κ of the dental X-ray test image g_(e)(n₁, n₂) to thedental X-ray reference image f_(e)(n₁, n₂) can be equal to 1 (one).Thereby, the positional displacement corrector 302 generates a dentalX-ray test image g_(eκθ)(n₁, n₂) which results from correcting therotation angle θ and the scaling ratio κ.

Subsequently, the positional displacement calculator 301 calculates aband-limited phase-only correlation function r_(fegeκθ) ^(K1K2)(n₁, n₂)between the contrast-corrected dental X-ray reference image f_(e)(n₁,n₂) and the dental X-ray test image g_(eκθ)(n₁, n₂) which results fromcorrecting the contrast, the rotation angle θ and the scaling ratio κ.In this respect, K₁/M₁=0.5 and K₂/M₂=0.5, for example. Thereafter, thepositional displacement calculator 301 detects the location of thecorrelation peak of the band-limited phase-only correlation functionr_(fegeκθ) ^(K1K2)(n₁, n₂), and thus calculates the amount of paralleltranslation of the dental X-ray test image g_(eκθ)(n₁, n₂) to the dentalX-ray reference image f_(e)(n₁, n₂).

The positional displacement corrector 302 corrects the dental X-ray testimage g_(eκθ)(n₁, n₂) which results from correcting the rotation angle θand the scaling ratio κ in order that the amount of parallel translationof the dental X-ray test image g_(eκθ)(n₁, n₂) to the dental X-rayreference image f_(e)(n₁, n₂) can be equal to 0 (zero), and thusgenerates a dental X-ray test image g′(n₁, n₂) and a dental X-rayreference image f′(n₁, n₂), shown in FIG. 8, which result fromcorrecting their respective rotation angles θ, scaling ratios κ andamounts of parallel translation.

The definer 303 shown in FIG. 1 defines any one of the dental X-ray testimage g′(n₁, n₂) and the dental X-ray reference image f′(n₁, n₂), whichare shown in FIG. 8, as a base image, and defines the other one of thetwo dental X-ray images as a corresponding image. The followingdescriptions will be provided with an assumption that the dental X-rayreference image f′(n₁, n₂) is defined as the base image whereas thedental X-ray test image g′(n₁, n₂) is defined as the correspondingimage. As shown in FIG. 9, the base point extractor 304 shown in FIG. 1extracts multiple base points from the base image f′(n₁, n₂) by usingthe Harris corner detector for detecting an abrupt change in the pixelvalues (or the brightness value). For the Harris corner detector,consider a matrix Q given by

$\begin{matrix}{Q = {\begin{matrix}\left\lbrack \frac{\partial I}{\partial n_{1}} \right\rbrack^{2} & {\frac{\partial I}{\partial n_{1}}\frac{\partial I}{\partial n_{2}}} \\{\frac{\partial I}{\partial n_{1}}\frac{\partial I}{\partial n_{2}}} & \left\lbrack \frac{\partial I}{\partial n_{2}} \right\rbrack^{2}\end{matrix}}} & {{Equation}\mspace{20mu} (5)}\end{matrix}$

where the pixel value at a point (n₁, n₂) is I(n₁, n₂).

When two eigenvalues of the matrix Q is large at a point in the baseimage f′(n₁, n₂), a corner exists there. A corner detector function E isgiven by

E=detQ−l(trQ)²   Equation (6)

where l is an adjustment parameter. A corner exists at a point to whichEquation (6) gives the locally largest values. In some cases, there arevery many points to which Equation (6) gives the locally largest valuesin the base image f′(n₁, n₂). With this taken into consideration, thebase point extractor 304 is designed to extract points, to which thelargest values not smaller than a predetermined threshold value aregive, as the base point. In this respect, a set of “B” base pointsextracted from the base image f′(n₁, n₂) is represented with U=(u₁*,u₂*, . . . , u_(B)*)^(T).

The corresponding point extractor 305, shown in FIG. 1, appliessub-pixel corresponding point search to a local image block around eachof the multiple base points, and thus finds corresponding points in thecorresponding image g′(n₁, n₂) as shown in FIG. 10. Specifically, thecorresponding point extractor 305, shown in FIG. 1, generates a roughbase image f′(n₁, n₂), which is rougher than the base image f′_(m)(n₁,n₂) by “m” degrees, in accordance with

$\begin{matrix}{{f_{m}^{\prime}\left( {n_{1},n_{2}} \right)} = {\frac{1}{4}{\sum\limits_{{p\; 1} = 0}^{m}{\sum\limits_{{p\; 2} = 0}^{m}{f_{m - 1}^{\prime}\left( {{{2n_{1}} + p_{1}},{{2n_{2}} + p_{2}}} \right)}}}}} & {{Equation}\mspace{20mu} (7)}\end{matrix}$

where m is an integer, and the base image f′(n₁, n₂) is represented asf′₀(n₁, n₂).

In addition, the corresponding point extractor 305 generates a roughcorresponding image g′_(m)(n₁, n₂), which is rougher than thecorresponding image g′(n₁, n₂) by “m” degrees, in accordance with

$\begin{matrix}{{g_{m}^{\prime}\left( {n_{1},n_{2}} \right)} = {\frac{1}{4}{\sum\limits_{{p\; 1} = 0}^{m}{\sum\limits_{{p\; 2} = 0}^{m}{g_{m - 1}^{\prime}\left( {{{2n_{1}} + p_{1}},{{2n_{2}} + p_{2}}} \right)}}}}} & {{Equation}\mspace{20mu} (8)}\end{matrix}$

where the corresponding image g′(n₁, n₂) is represented as g′₀(n₁, n₂).

In Equations (7) and (8), the largest value of m is set at 2, forexample. Subsequently, the corresponding point extractor 305 calculatesa rough base point u_(m)=(u_(m 1), u_(m 2)) in the rough base imagef′_(m)(n₁, n₂), which corresponds to a base point u=u₀=(u₁, u₂) in thebase image f′₀(n₁, n₂), in accordance with

u _(m)=[(1/2)u _(m−1)]=([(1/2)u _(m−·1)],[(1/2)u _(m−1·2)])   Equation(9)

where [x] means that the fractional part of x is rounded off.

Furthermore, when “m” takes on its maximum, the corresponding pointextractor 305 assumes that the coordinate of the rough correspondingpoint v_(m) in the rough corresponding image g′_(m)(n₁, n₂) is equal tothe coordinate of the rough base point u_(m) in the rough base imagef′_(m)(n₁, n₂).

Subsequently, the corresponding point extractor 305 extracts a sub roughbase image f′_(m−1sub)(n₁, n₂) with the center on the rough base pointu_(m−1) from the rough base image f′_(m−1)(n₁, n₂). In addition, thecorresponding point extractor 305 extracts a sub rough correspondingimage g′_(m−1sub)(n₁, n₂) around a coordinate 2v_(m) from the roughcorresponding image g′_(m−1)(n₁, n₂). The size of each of the sub roughbase image f′_(m−1sub)(n₁, n₂) and the sub rough corresponding imageg′_(m−1sub)(n₁, n₂) is 32×32 pixels, for example.

Thereafter, the corresponding point extractor 305 calculates apositional displacement between the sub rough base image f′_(m−1sub)(n₁,n₂) and the sub rough corresponding image g′_(m−1sub)(n₁, n₂) by usingthe window function and the phase-only correlation function given by

$\begin{matrix}{{w\left( {n_{1},n_{2}} \right)} = {\frac{1 + {\cos \left( {\pi \; {n_{1}/M_{1}}} \right)}}{2}{\frac{1 + {\cos \left( {\pi \; {n_{2}/M_{2}}} \right)}}{2}.}}} & {{Equation}\mspace{20mu} (10)}\end{matrix}$

Note that the use of the window function solves the discontinuity of theimage ends.

When the calculated positional displacement vector is denoted byδ_(m−1), a rough corresponding point v_(m−1) in the sub roughcorresponding image g′_(m−1sub)(n₁, n₂), which corresponds to the roughbase point u_(m−1) in the sub rough base image f′_(m−1sub)(n₁, n₂), isgiven by

ν_(m−1)=2ν_(m)+δ_(m−1).   Equation (11)

The corresponding point extractor 305 calculates a rough correspondingpoint v₀ by sequentially decreasing “m” one-by-one from its maximum.When, for example, the maximum of “m” is equal to 2, the correspondingpoint extractor 305 calculates a rough base point u₂=(u₂ ₁, u₂ ₂) in arough base image f′₂(n₁, n₂) by using Equation (9), and thus assumesthat a rough corresponding point v₂ in a rough corresponding imageg′₂(n₁, n₂) is equal to u₂. Subsequently, the corresponding pointextractor 305 extracts a sub rough base image f′_(1sub)(n₁, n₂) with thecenter thereof on a rough base point u₁ from a rough base image f′₁(n₁,n₂), and extracts a sub rough corresponding image g′_(1sub)(n₁, n₂) withthe center thereof on a point 2v₂ from a rough corresponding imageg′₁(n₁, n₂). In addition, the corresponding point extractor 305calculates a positional displacement δ₁ between the sub rough base imagef′_(1sub)(n₁, n₂) and the sub rough corresponding image g′_(1sub)(n₁,n₂) by using the window function and the phase-only correlationfunction, and calculates a rough corresponding point v₁ by usingEquation (11).

Thereafter, the corresponding point extractor 305 extracts a sub baseimage f′_(0sub)(n₁, n₂) around a base point u₀ from a base image f′₀(n₁,n₂), and extracts a sub corresponding image g′_(0sub)(n₁, n₂) around apoint 2v₁ from a corresponding image g′₀(n₁, n₂). Furthermore, thecorresponding point extractor 305 calculates a positional displacementδ₀ between the sub base image f′_(0sub)(n₁, n₂) and the subcorresponding image g′_(0sub)(n₁, n₂) by using the window function andthe phase-only correlation function, and calculates a roughcorresponding point v₀=v=(v₁, v₂) by using Equation (11). The calculatedrough corresponding point v₀=v=(v₁, v₂) in the corresponding imageg′₀(n₁, n₂) includes no pixel-level error.

Thereafter, the corresponding point extractor 305 extracts the sub baseimage f′_(0sub)(n₁, n₂) around a base point u=u₀=(u₁, u₂) from the baseimage f′₀(n₁, n₂). In addition, the corresponding point extractor 305extracts the sub corresponding image g′_(0sub)(n₁, n₂) around the roughcorresponding point v₀=v=(v₁, v₂) calculated from the correspondingimage g′₀(n₁, n₂) at a pixel-level. In this respect, n₁=−M_(1sub), . . .,M_(1sub)(M_(1sub)>0), n₂=−M_(2sub), . . . ,M_(2sub)(M_(2sub)>0),N_(1sub)=2_(M) _(1sub)+1, and N_(2sub)=2M_(2sub)+1.

Furthermore, the corresponding point extractor 305 defines an initialwindow function w_(f)(n₁, n₂) for the sub base image f′_(0sub)(n₁, n₂)and an initial window function w_(g)(n₁, n₂) for the sub correspondingimage g′_(0sub)(n₁, n₂). A window function w(x₁, x₂) corresponding toall of the real variables x₁ and x₂ is given by

$\begin{matrix}{{w\left( {x_{1},x_{2}} \right)} = {\frac{1 + {\cos \left( {\pi \; {x_{1}/M_{1\; {sub}}}} \right)}}{2}\frac{1 + {\cos \left( {\pi \; {x_{2}/M_{2{sub}}}} \right)}}{2}}} & {{Equation}\mspace{20mu} (12)}\end{matrix}$

where w_(f)(n₁, n₂)=w_(g)(n₁, n₂)=w(x₁, x₂)|_(x1=n1, x2=n2).

Instead of the Hanning window function given by Equation (12), theHamming, Gaussian or Kaiser window function can be used.

Thereafter, the corresponding point extractor 305 calculates apositional displacement δ=(δ₁, δ₂)between the sub base imagef′_(0sub)(n₁, n₂) and the sub corresponding image g′_(0sub)(n₁, n₂) byusing the initial window function and the phase-only correlationfunction. In this respect, when a fine positional displacement δ=(δ₁,δ₂) between the sub base image f′_(0sub)(n₁, n₂) and the subcorresponding image g′_(0sub)(n₁, n₂) is represented with δ=(δ₁, δ₂),the phase-only correlation function of the sub base image f′_(0sub)(n₁,n₂) and the sub corresponding image g′_(0sub)(n₁, n₂) is modeled by

$\begin{matrix}{{r\left( {n_{1},n_{2}} \right)} \approx {\frac{\alpha}{N_{1{sub}}N_{2{sub}}}\frac{\sin \left\{ {\pi \left( {n_{1} + \delta_{1}} \right)} \right\}}{\sin \left\{ {\frac{\pi}{N_{1{sub}}}\left( {n_{1} + \delta_{1}} \right)} \right\}}\frac{\sin \left\{ {\pi \left( {n_{2} + \delta_{2}} \right)} \right\}}{\sin \left\{ {\frac{\pi}{N_{2{sub}}}\left( {n_{2} + \delta_{2}} \right)} \right\}}}} & {{Equation}\mspace{20mu} (13)}\end{matrix}$

where α=1. α is a parameter introduced to represent the height of thecorrelation peak. α is actually not larger than 1 (one), because thevalue of α decreases when uncorrelated noise is added to the images. Byfitting the correlation peak model given by Equation (13) to the numericdata of the actually computed phase-only correlation function, it ispossible to estimate the peak location existing between the pixels, andaccordingly to estimate the parameters α, δ₁ and δ₂. The estimated δ₁and δ₂ are smaller than 1 (one), and are at a sub-pixel level.Furthermore, the corresponding point extractor 305 updates the windowfunctions w_(f)(n₁, n₂) and w_(g)(n₁, n₂), with the positionaldisplacement δ=(δ₁, δ₂) calculated at a sub-pixel level, by use of

$\begin{matrix}{{w_{f}\left( {n_{1},n_{2}} \right)} = \left. {w\left\lbrack {{x_{1} + \frac{\delta_{1}}{2}},{x_{2} + \frac{\delta_{2}}{2}}} \right\rbrack} \right|_{{{x\; 1} = {n\; 1}},{{x\; 2} = {n\; 2}}}} & {{Equation}\mspace{20mu} (14)} \\{and} & \; \\{{w_{g}\left( {n_{1},n_{2}} \right)} = \left. {w\left\lbrack {{x_{1} - \frac{\delta_{1}}{2}},{x_{2} - \frac{\delta_{2}}{2}}} \right\rbrack} \middle| {}_{{{x\; 1} = {n\; 1}},{{x\; 2} = {n\; 2}}}. \right.} & {{Equation}\mspace{20mu} (15)}\end{matrix}$

Subsequently, the corresponding point extractor 305 repeats thecalculation of the positional displacement δ=(δ₁, δ₂) by using theupdated window functions and the phase-only correlation function. Thenumber of times the corresponding point extractor 305 repeats thecalculation is five, for example. Each time the calculation of thepositional displacement δ=(δ₁, δ₂) is repeated, the sub-pixel-levelerror becomes smaller. Thereafter, the corresponding point extractor 305adds the positional displacement δ=(δ₁, δ₂) to the rough correspondingpoint v₀=v=(v₁, v₂) as shown by

ν_(0F)=ν₀+δ  Equation (16)

and thus calculates a corresponding point v_(0f) in the correspondingimage g′(n₁, n₂).

In this respect, a set of “B” corresponding points found in thecorresponding image g′(n₁, n₂) is represented with V=(v₁*, v₂*, . . . ,v_(B)*)^(T). Note that the corresponding point extractor 305 is designedto find a corresponding point for each of the multiple base points.However, if the correlation peak value of the phase-only correlationfunction is smaller than the threshold value, the corresponding pointsthus found are eliminated as false corresponding points. The thresholdvalue is 0.4, for example.

The correspondence calculator 306, shown in FIG. 1, holds the thin platespline (TPS) function given by

v=TPS(u)=d+A·u+C ^(T) ·s(u)   Equation (17)

where “u” denotes a pixel in the base image f′(n₁, n₂), “v” denotes apixel in the corresponding image g′(n₁, n₂), “d” denotes a 2×1translation vector, “A” denotes a 2×2 affine translation matrix, and “C”denotes a “B×2” coefficient matrix. In addition, s(u)=(σ(u−u₁*),σ(u−u₂*), . . . , σ(u−u_(B)*))^(T). σ(u) is represented with

$\begin{matrix}{{\sigma (u)} = \left\{ \begin{matrix}{{u}^{2}{\log \left( {u} \right)}} & {{u} > 0} \\0 & {{u} = 0.}\end{matrix} \right.} & {{Equation}\mspace{20mu} (18)}\end{matrix}$

Equation (17) has 2B+6 parameters. The correspondence calculator 306substitutes a set U=(u₁*, u₂*, . . . , u_(B)*)^(T) of “B” base pointsand a set V=(v₁*, v₂*, . . . , v_(B)*)^(T) of “B” corresponding pointsinto Equation (19), and finds 2B+6 parameters by solving Equation (19)using least square method.

$\begin{matrix}{{\begin{bmatrix}H & l_{B} & U \\l_{B}^{T} & 0 & 0 \\U^{T} & 0 & 0\end{bmatrix}\begin{bmatrix}C \\d^{T} \\A^{T}\end{bmatrix}} = \begin{bmatrix}V \\0 \\0\end{bmatrix}} & {{Equation}\mspace{20mu} (19)}\end{matrix}$

where 1_(B) denotes a vector consisting of “l” with a length “B”, “H”denotes a B×B matrix with an element h_(ij)=σ(u_(i)*−u_(j)*).Subsequently, the correspondence calculator 306 solves Equation (19) byusing the method of least square, and thereby finds 2B+6 parameters.

The correspondence calculator 306 substitutes the 2B+6 parameters thusfound into Equation (17), and calculates the thin plate spline functionrepresenting correspondence between the multiple base points included inthe base image and the multiple corresponding points included in thecorresponding image. The nonlinear distortion corrector 307 transformsthe corresponding image g′(n₁, n₂), as shown by grids in FIG. 11, byusing the thin plate spline function representing the correspondencebetween the multiple base points included in the base image and themultiple corresponding points included in the corresponding image.Thereafter, the nonlinear distortion corrector 307, shown in FIG. 1,converts all of the coordinates constituting the corresponding imageg′(n₁, n₂), based on the Equation (17) into which the parameters aresubstituted, and thus generates a corresponding image g″(n₁, n₂) whichresults from correcting its nonlinear distortion from the base imagef′(n₁, n₂), as shown in FIG. 12.

The CPU 300, shown in FIG. 1, further includes a common area extractor402. An area and coordinates which is included in any one of the baseimage f′(n₁, n₂) and the corresponding image g″(n₁, n₂) but not in theother one of the two images represent noise components uncorrelated tothe phase-only correlation function. As shown in FIG. 13, the commonarea extractor 402 extracts a common area between the base image f′(n₁,n₂) and the corresponding image g″(n₁, n₂) which results from correctingits nonlinear distortion from the base image f′(n₁, n₂), and thusgenerates a common area f′″(n₁, n₂) extracted from the base image f′(n₁,n₂) and a common area g′″(n₁, n₂) extracted from the corresponding imageg″(n₁, n₂). The common area extractor 402 extracts the common areas byusing a pixel projection in an n₁ direction and a pixel projection in ann₂ direction, respectively. The size of the common area f′″(n₁, n₂)extracted from the base image f′(n₁, n₂) is equal to the size of thecommon area g′″(n₁, n₂) extracted from the corresponding image g″(n₁,n₂).

The similarity calculator 308, shown in FIG. 1, calculates theband-limited phase-only correlation function of the common area f′″(n₁,n₂) and the common area g′″(n₁, n₂) with a condition of, for example,K₁/M₁=K₂/M₂=0.1. Furthermore, the similarity calculator 308 finds thelargest peak value of the band-limited phase-only correlation functionas a comparison score for indicating a similarity between the dentalX-ray test image g(n₁, n₂) and the dental X-ray reference image f(n₁,n₂). The CPU 300 further includes an identifier 403. If a comparisonscore indicating a similarity between the dental X-ray test image g(n₁,n₂) of a photographed person who is not identified and the dental X-rayreference image f(n₁, n₂) of a photographed person which is identifiedis higher than a threshold value, the identifier 403 identifies theperson of whom the dental X-ray test image g(n₁, n₂) is taken as theperson of whom the dental X-ray reference image f(n₁, n₂) is taken.

A data memory 201 is connected to the CPU 300. The data memory 201includes a reference image memory module 202, a personal informationmemory module 203 and an identification result memory module 204. Thedental X-ray reference image f(n₁, n₂) is stored in the reference imagememory module 202. Information on a person whom the dental X-rayreference image f(n₁, n₂) has been taken of is stored in the personalinformation memory module 203. A result of identification carried out bythe identifier 403 is stored in the identification result memory module204.

An input device 312, an output device 313, a program memory 330 and atemporary memory 331 are further connected to the CPU 300. Examples ofwhat is usable as the input device 312 include a keyboard and a pointingdevice such as a mouse. Examples of what is usable as the output device313 include image displaying devices such as a liquid crystal displayand a monitor, and a printer. An operating system and the like forcontrolling the CPU 300 are stored in the program memory 330. Results ofarithmetic performed by the CPU 300 are sequentially stored in thetemporary memory 331. Examples of what are usable as the program memory330 and the temporary memory 331 include storage media, such as asemiconductor memory, a magnetic disc, an optical disc, amagneto-optical disc and a magnetic tape, in which a program is stored.

Next, descriptions will be provided for a method for comparing dentalX-ray images according to the embodiment by using a flowchart shown inFIG. 14. It should be noted that results of arithmetic performed by theCPU 300 shown in FIG. 1 are sequentially stored in the temporary memory331.

In step S101, the contrast corrector 401 reads the dental X-ray testimage g(n₁, n₂), shown in FIG. 6, through the image acquiring device200. In addition, the contrast corrector 401, shown in FIG. 1, reads thedental X-ray reference image f(n₁, n₂), shown in FIG. 6, which has beensaved in the reference image memory module 202. Subsequently, thecontrast corrector 401 generates a contrast-corrected dental X-ray testimage g_(e)(n₁, n₂) and a contrast-corrected dental X-ray referenceimage f_(e)(n₁, n₂), which are shown in FIG. 7. Thereafter, the contrastcorrector 401, shown in FIG. 1, transmits the contrast-corrected dentalX-ray test image g_(e)(n₁, n₂) and the contrast-corrected dental X-rayreference image f_(e)(n₁, n₂) to the positional displacement calculator301.

In step S102, the positional displacement calculator 301 calculates thetwo-dimensional discrete Fourier transforms G_(e)(k₁, k₂) and F_(e)(k₁,k₂) respectively of the dental X-ray test image g_(e)(n₁, n₂) and thedental X-ray reference image f_(e)(n₁, n₂) in accordance with Equations(1) and (2). Subsequently, the positional displacement calculator 301calculates G_(LP)(k₁′, k₂′) resulting from the log-polar transform ofthe square root of the amplitude spectrum of the dental X-ray test image(|G_(e)(k₁, k₂)|)^(1/2) and F_(LP)(k₁′, k₂′) resulting from thelog-polar transform of the square root of the amplitude spectrum of thedental X-ray reference image (|F_(e)(k₁, k₂)|)^(1/2). Thereafter, thepositional displacement calculator 301 calculates the normalized crosspower spectrum R_(FlpGlp)(k₁′, k₂′) in accordance with Equation (3), andthereafter calculates the band-limited phase-only correlation functionr_(FlpGlp) ^(K1K2)(n₁, n₂) in accordance with Equation (4).

In step S103, the positional displacement calculator 301 calculates therotation angle θ of the dental X-ray test image g_(e)(n₁, n₂) to thedental X-ray reference image f_(e)(n₁, n₂) and the scaling ratio κ ofthe dental X-ray test image g_(e)(n₁, n₂) to the dental X-ray referenceimage f_(e)(n₁, n₂), based on the location of the correlation peak ofthe band-limited phase-only correlation function r_(FlpGlp) ^(K1K2)(n₁,n₂), and transmits the rotation angle θ and the scaling ratio κ to thepositional displacement corrector 302. The positional displacementcorrector 302 corrects the dental X-ray test image g_(e)(n₁, n₂) inorder that the rotation angle θ can become equal to 0 (zero) and inorder that the scaling ratio κ can become equal to 1 (one), and thusgenerates the dental X-ray test image g_(eκθ)(n₁, n₂) which results fromcorrecting the rotation angle θ and the scaling ratio κ. The positionaldisplacement corrector 302 transmits the corrected dental X-ray testimage g_(eκθ)(n₁, n₂) to the positional displacement calculator 301.

In step S104, the positional displacement calculator 301 calculates theband-limited phase-only correlation function r_(fegeκθ) ^(K1K2)(n₁, n₂)between the contrast-corrected dental X-ray reference image f_(e)(n₁,n₂) and the dental X-ray test image g_(eκθ)(n₁, n₂) which results fromcorrecting the contrast, the rotation angle θ and the scaling ratio κ.Subsequently, the positional displacement calculator 301 calculates theamount of parallel translation of the dental X-ray test imageg_(eκθ)(n₁, n₂) which results from correcting the rotation angle θ andthe scaling ratio κ of the dental X-ray test image to the dental X-rayreference image f_(e)(n₁, n₂), based on the location of the correlationpeak. Thereafter, the positional displacement calculator 301 transmitsthe amount of parallel translation thus calculated to the positionaldisplacement corrector 302.

In step S105, the positional displacement corrector 302 generates thedental X-ray test image g′(n₁, n₂) and the dental X-ray reference imagef′(n₁, n₂), shown in FIG. 8, to which a correction is applied in orderthat their respective amount of parallel translation can become equal to0 (zero). The positional displacement corrector 302, shown in FIG. 1,transmits the generated dental X-ray test image g′(n₁, n₂) and thegenerated dental X-ray reference image f′(n₁, n₂) to the definer 303. Instep S106, the definer 303 defines the dental X-ray reference imagef′(n₁, n₂) as the base image, and defines the dental X-ray test imageg′(n₁, n₂) as the corresponding image. It should be noted that thedefiner 303 may define the dental X-ray reference image f′(n₁, n₂) asthe corresponding image, and the dental X-ray test image g′(n₁, n₂) asthe base image. Subsequently, the definer 303 transmits the base imagef′(n₁, n₂) to the base point extractor 304, and transmits thecorresponding image g′(n₁, n₂) to the corresponding point extractor 305.In step S107, as shown in FIG. 9, the base point extractor 304 extractsthe set U=(u₁*, u₂*, . . . , u_(B)*)^(T) of “B” base points from thebase image f′(n₁, n₂) by using the Harris corner detection usingEquations (5) and (6) or the like. Thereafter, the base point extractor304, shown in FIG. 1, transmits the base image f′(n₁, n₂) and the setU=(u₁*, u₂*, . . . , u_(B)*)^(T) of the base points to the correspondingpoint extractor 305.

In step S108, the corresponding point extractor 305 calculates themultiple rough corresponding points which respectively correspond to themultiple base points extracted from the base image f′(n₁, n₂) by usingEquations (7) to (11). Subsequently, the corresponding point extractor305 calculates the positional displacement δ=(δ₁, δ₂) between the subbase image f′_(0sub)(n₁, n₂) around each base point and the subcorresponding image g′_(0sub)(n₁, n₂) around each rough correspondingpoint by using the initial window function given by Equation (12) andthe phase-only correlation function. Thereafter, the corresponding pointextractor 305 updates the window functions w_(f)(n₁, n₂) and w_(g)(n₁,n₂) as shown by Equations (14) and (15), and thus repeats thecalculation of the positional displacement δ=(δ₁, δ₂), and accordinglydetermines the corresponding point corresponding to the base point. Thecorresponding point extractor 305 determines the corresponding point foreach of the base points, and thus determines the set V=(v₁*, v₂*, . . ., v_(B)*)^(T) of the B corresponding points, which corresponds to theset U=(u₁*, u₂*, . . . , u_(B)*)^(T) of the base points, in thecorresponding image g′(n₁, n₂), as shown in FIG. 10. After that, thecorresponding point extractor 305, shown in FIG. 1, transmits the setU=(u₁*, u₂*, . . . , u_(B)*)^(T) of the base points and the set V=(v₁*,v₂*, . . . , v_(B)*)^(T) of the corresponding points to thecorrespondence calculator 306.

In step S109, the correspondence calculator 306 finds the 2B+6parameters of the thin-plate spline function given by Equation (17) byusing the set U=(u₁*, u₂*, . . . , u_(B)*)^(T) of the base points andthe set V=(v₁*, v₂*, . . . , v_(B)*)^(T) of the corresponding points inaccordance with Equation (19). Subsequently, the correspondencecalculator 306 substitutes the 2B+6 parameters thus found into Equation(17). Thereafter, the correspondence calculator 306 transmits thethin-plate spline function into which the 2B+6 parameters aresubstituted to the nonlinear distortion corrector 307. In step S110, thenonlinear distortion corrector 307 converts all of the coordinatesconstituting the corresponding image g′(n₁, n₂), based on the Equation(17) into which the 2B+6 parameters are substituted, and thus generatesthe corresponding image g″(n₁, n₂) which results from correcting thenonlinear distortion of the corresponding image g′(n₁, n₂) to the baseimage f′(n₁, n₂), as shown in FIG. 12. After that, the nonlineardistortion corrector 307, shown in FIG. 1, transmits the correspondingimage g″(n₁, n₂), which results from correcting the nonlinear distortionof the corresponding image g′(n₁, n₂), to the common area extractor 402.In step S111, as shown in FIG. 13, the common area extractor 402generates the common area f′″(n₁, n₂) extracted from the base imagef′(n₁, n₂) and the common area g′″(n₁, n₂) extracted from thecorresponding image g″(n₁, n₂). The nonlinear distortion corrector 307,shown in FIG. 1, transmits the common area f′″(n₁, n₂) and the commonarea g′″(n₁, n₂) to the similarity calculator 308.

In step S112, the similarity calculator 308 finds the largest peak valueof the band-limited phase-only correlation function of the common areaf′″(n₁, n₂) and the common area g′″(n₁, n₂) as the comparison score forindicating the similarity between the dental X-ray test image g(n₁, n₂)and the dental X-ray reference image f(n₁, n₂). Subsequently, thesimilarity calculator 308 transmits the comparison score to theidentifier 403. Instep S113, the identifier 403 reads from the personalinformation memory module 203 the information on the person of whom thedental X-ray reference image f(n₁, n₂) has been taken. Thereafter, theidentifier 403 identifies the person of whom the dental X-ray test imageg(n₁, n₂) has been taken as the person of whom the dental X-rayreference image f(f₁, f₂) has been taken, when the comparison score ishigher than the threshold value. After that, the identifier 403 storesin the identification result memory module 204 the information that theperson of whom the dental X-ray test image g(n₁, n₂) has been taken isidentified as the person of whom the dental X-ray reference image f(f₁,f₂) has been taken. With this, the identifier 402 ends the method ofcomparing dental X-ray images according to the embodiment.

The above-described system and method for comparing dental X-ray images,according to the embodiment, make it possible to compare the dentalX-ray test image g(n₁, n₂) and the dental X-ray reference image f(f₁,f₂) with a higher accuracy, even in a case where a nonlinear distortionexists between the dental X-ray test image g(n₁, n₂) and the dentalX-ray reference image f(f₁, f₂).

EXAMPLE

Similarities were evaluated between dental radiographs taken before therespective dental treatments and dental radiographs taken after therespective dental treatments by using a database which stores dentalradiographs taken before and after dental treatments. For this example,the database was constructed as follows. A dental radiograph (with367×485 pixels) was taken of each of 250 patients before their dentaltreatments. A week or more later than their dental treatments, anotherdental radiograph (with 367×485 pixels) was taken of each of the same250 patients. The total of 500 radiographs (250 patients×2 radiographs)were stored in the database. Note that, because each dental radiographhas a margin in each side, an image with 307×425 pixels obtained byexcluding 30 pixels in each side was used for each dental radiograph.Each of FIGS. 15, 16 and 17 shows examples respectively of a dentalradiograph taken before a dental treatment and another dental radiographtaken after the dental treatment, which were stored in the database.

The method for comparing dental X-ray images according to the embodimentand two methods each of comparing dental X-ray images according tocomparative examples were used for the similarity evaluation. In thecase of the comparing method according to the first comparative example,a correction was applied to the dental radiographs taken before andafter the dental treatments in terms of a rotation angle θ, a scalingratio κ and the amount of parallel translation only, but not in terms ofa nonlinear distortion. In the case of the comparing method according tothe second comparative example, a correction was applied to the dentalradiographs taken before and after the dental treatments in terms of arotation angle θ, a scaling ratio κ and the amount of paralleltranslation, and subsequently another correction was applied to the samedental radiographs in terms of a nonlinear distortion by using theprojection transform.

For quantitative evaluation of the positioning accuracy, verificationwas made on a 1-to-n basis. In addition, the dental radiographs takenafter the respective treatments were defined as “dental X-ray testimages,” and the dental radiographs taken before the respectivetreatment which were stored in the database were defined as “dentalX-ray reference images.” One “dental X-ray test image” was compared withthe 250 “dental X-ray reference images.” Thereby, the comparison scorewas calculated for each paired dental X-ray images. This process wasapplied to all of the “dental X-ray test images.”

Subsequently, the rank of each comparison score was plotted on thehorizontal axis. At each rank, a verification probability representing aprobability that the “dental X-ray test image” and the “dental X-rayreference image” may be taken of the same person was plotted on thevertical axis. Thereby, a cumulative match curve (CMC) shown in FIG. 18was obtained as follows. The cumulative match curve is a scheme forevaluating the performance of the 1-to-n verification. A higherverification probability with a higher comparison score means that theimages were correctly matched.

When the method for comparing dental X-ray images according to theembodiment was used, the probability that a pair of the “dental X-raytest image” and the “dental X-ray reference image” both taken of thesame person recorded the first highest comparison score was 75.6%(=189/250). In contrast, when the method according to the firstcomparative example was used, the probability that a pair of the “dentalX-ray test image” and the “dental X-ray reference image” both taken ofthe same person recorded the first highest comparison score was 59.2%(=148/250). When the method according to the second comparative examplewas used, the probability that a pair of the “dental X-ray test image”and the “dental X-ray reference image” both taken of the same personrecorded the first highest comparison score was 66.8% (=167/250).

Furthermore, when the method for comparing dental X-ray images accordingto the embodiment was used, all of the pairs each consisting of the“dental X-ray test image” and the “dental X-ray reference image” bothtaken of the same person were able to be included in the pairs whichrecorded the first to 14th highest comparison scores for each dentalX-ray test image. By contrast, when the methods according to the firstand second comparative examples were used, not all of the pairs eachconsisting of the “dental X-ray test image” and the “dental X-rayreference image” both taken of the same person were able to be includedin the pairs which recorded the first to 20th highest comparison scoresfor each dental X-ray test image.

Moreover, as shown in FIG. 19, generally, the comparison score of eachpaired “dental X-ray test image” and “dental X-ray reference image” bothtaken of the same person was higher when the method of comparing dentalX-ray images according to the embodiment was used than when the methodsaccording to the first and second comparative examples were used.

FIGS. 20 to 23 each show the “dental X-ray reference image”, the “dentalX-ray test image,” a “dental X-ray test image” corrected by using themethod according to the first comparative example, a “dental X-ray testimage” corrected by using the method according to the second comparativeexample, a “dental X-ray test image” and “dental X-ray reference image”corrected by using the method of comparing dental X-ray images accordingto the embodiment, as well as a differential image of the “dental X-raytest image” corrected by using the method of comparing dental X-rayimages according to the embodiment. As shown in FIGS. 20 to 23, the useof the method of comparing dental X-ray images according to theembodiment made the comparison score higher than the use of any othermethod.

Other Embodiments

The foregoing descriptions have been provided for the present inventionciting the embodiment. However, the descriptions and drawingsconstituting this disclosure should not be understood as limiting thepresent invention. From this disclosure, various alternativeembodiments, examples and operational techniques will be clear to thoseskilled in the art. For example, it is possible to identify theunidentified dead body with a higher accuracy by taking a dental X-raytest image g(n₁, n₂) of the dead body, and thus by comparing the dentalX-ray test image g(n₁, n₂) with dental X-ray reference images f(n₁, n₂)taken of identified persons. It should be understood that, in thismanner, the present invention includes various embodiments and the likewhich have not been described herein. For this reason, the presentinvention is limited only by appropriate matters defining the inventionin the scope of claims.

1. A system for comparing dental X-ray images comprising: a positionaldisplacement calculator configured to calculate a positionaldisplacement between a dental X-ray test image and a dental X-rayreference image, by using phase-only correlation; a positionaldisplacement corrector configured to correct the positionaldisplacement; a base point extractor configured to define, as a baseimage, any one of the dental X-ray test image and the dental X-rayreference image between which the positional displacement is corrected,to define, as a corresponding image, the other one of the two dentalimages, and to extract a plurality of base points from the base image; acorresponding point extractor configured to extract a plurality ofcorresponding points, which correspond to the plurality of base points,from the corresponding image; a correspondence calculator configured tocalculate correspondence between the plurality of base points and theplurality of corresponding points nonlinearly distorted from theplurality of base points; a nonlinear distortion corrector configured tocorrect a nonlinear distortion between the base image and thecorresponding image, based on the correspondence; and a similaritycalculator configured to find, by using phase-only correlation, asimilarity between the base image and the corresponding image betweenwhich the nonlinear distortion is corrected.
 2. The system of claim 1,wherein the correspondence is given by a thin plate spline function. 3.The system of claim 1, wherein a pixel value changes beyond a thresholdvalue at each of the plurality of base points.
 4. The system of claim 1,wherein the phase-only correlation is a band-limited phase-onlycorrelation function.
 5. The system of claim 1, further comprising acontrast corrector configured to correct a contrast of the dental X-raytest image.
 6. The system of claim 1, further comprising a contrastcorrector configured to correct a contrast of the dental X-ray referenceimage.
 7. The system of claim 1, further comprising a common areaextractor configured to extract a common area between the base image andthe corresponding image between which the nonlinear distortion iscorrected.
 8. The system of claim 1, wherein a person whom the dentalX-ray test image has been taken of is unidentified, and a person whomthe dental X-ray reference image has been taken of is identified.
 9. Thesystem of claim 8, further comprising a personal information memorymodule configured to store information on the person whom the dentalX-ray reference image has been taken of.
 10. The system of claim 9,further comprising an identifier configured to identify the person whomthe dental X-ray test image has been taken of as the person whom thedental X-ray reference image has been taken of when the base image andthe corresponding image are similar to each other.
 11. The system ofclaim 1, wherein the similarity calculator uses, as a comparison score,a maximum peak value of the phase-only correlation function used for thephase-only correlation.
 12. A method for comparing dental X-ray imagesincluding: calculating a positional displacement between a dental X-raytest image and a dental X-ray reference image, by using phase-onlycorrelation; correcting the positional displacement; defining, as a baseimage, any one of the dental X-ray test image and the dental X-rayreference image between which the nonlinear distortion is corrected, anddefining, as a corresponding image, the other one of the two dentalimages, as well as extracting a plurality of base points from the baseimage; extracting a plurality of corresponding points, which correspondto the plurality of base points, from the corresponding image;calculating correspondence between the plurality of base points and theplurality of corresponding points nonlinearly distorted from coordinatesof the plurality of base points; correcting a nonlinear distortionbetween the base image and the corresponding image, based on thecorrespondence; and finding, by using phase-only correlation, asimilarity between the base image and the corresponding image betweenwhich the nonlinear distortion is corrected.
 13. The method of claim 12,wherein the correspondence is given by a thin plate spline function. 14.The method of claim 12, wherein a pixel value changes beyond a thresholdvalue at each of the plurality of base points.
 15. The method of claim12, wherein the phase-only correlation is a band-limited phase-onlycorrelation function.
 16. The method of claim 12 further includingcorrecting a contrast of the dental X-ray test image.
 17. The method ofclaim 12 further including correcting a contrast of the dental X-rayreference image.
 18. The method of claim 12 further including extractinga common area between the base image and the corresponding image betweenwhich the nonlinear distortion is corrected.
 19. The method of claim 12,wherein a person whom the dental X-ray test image has been taken of isunidentified, and a person whom the dental X-ray reference image hasbeen taken of is identified.
 20. The method of claim 19 furtherincluding identifying the person whom the dental X-ray test image hasbeen taken of as the person whom the dental X-ray reference image hasbeen taken of when the base image and the corresponding image aresimilar to each other.
 21. The method of claims 12 to 20, wherein, inthe similarity finding step, a maximum peak value of the phase-onlycorrelation function used for the phase-only correlation is used as acomparison score.